arXiv:cs/0512072vl [cs.SC] 18 Dec 2005 


Computations with one and two real algebraic numbers 


Ioannis Z. Emir is and Elias P. Tsigaridas 
Department of Informatics and Telecommunications 
National University of Athens, HELLAS 
{emiris,et}@di.uoa.gr 


Abstract 

We present algorithmic and complexity results concerning computations with one and 
two real algebraic numbers, as well as real solving of univariate polynomials and bivariate 
polynomial systems with integer coefficients using Sturm-Habicht sequences. 

Our main results, in the univariate case, concern the problems of real root isolation fTh. 1191) 
and simultaneous inequalities (Cor. 1261 and in the bivariate, the problems of system real 
solving ITh. 1421 . sign evaluation ITh. 1371 and simultaneous inequalities (Cor. 031 ■ 


1 Introduction 

In what follows O b means bit complexity and the 0 b- notation means that we are ignoring loga¬ 
rithmic factors. 

For a polynomial P £ Z[X,Y], deg(P) denotes its total degree while deg x (P) (resp. deg Y (P)) 
denotes its degree if we consider it as a univariate polynomial with respect to X (resp. Y). 

By £ (P) we denote an upper bound on the bit size of the coefficients of P (including a bit for 
the sign) i.e £ (P) = Llg(max |at|)J + 2, where at are the coefficients of P. In particular £ (a) is 
the bit size of a, if it is a non zero integer, or the maximum bit size of the numerator and the 
denominator if it is a rational number. 

Let M (t) denote the bit complexity of multiplying two integers of bit size at most x and 
M (d, x) denote the bit complexity of multiplying two univariate polynomials of degrees bounded 
by d and coefficient bit size at most x. Using fast multiplication algorithms, i.e FFT umm, 
the complexities of these operations are 

M (x) = D(xlgxlglgx) 

M(d,x) = 0 B (dxlg(dx) lglg(dx)) ^ 

Hence M (x) = 0 B (xlg c ' x) and M (d,x) = 0 B (dxlg C2 (dx)) for suitable constants c^ ,C 2 - 


2 Sturm-Habicht sequences 

Let A = H B= o Q kX k ,B = )Lk=o bkX k £ Z[X] where deg(A) = p > q = deg(B) and £ (A) = 
£ (B) = x. 

We denote by rem(A,B) and quo (A,B) the remainder and the quotient, respectively, of the 
division of A and B. In general both may have rational coefficients. Following J2S1 we give the 
following definitions: 

Definitions 1 The signed polynomial remainder sequence of A and B. sPRS (A, B) is the poly¬ 
nomial sequence 

(R 0 = A, R] = B, R 2 = — rem (A, B) ,..., R k = — rem (Rk-2, Rk-i)) 
where rem (Rk-i, Rk) = 0. 
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The quotient sequence of A and B is the polynomial sequence 

(Qo = quo (R 0 , Ri ) , Qi =quo(Ri,R 2 ),...,Qk-i = quo (R k -i , R k )) 

We also need the definition of the quotient boot, which is the polynomial sequence 


(Qo, Qi, ■ • •, Qk-i, Rk) 


There is a huge bibliography on signed polynomial remainder sequences (c.f [annum and 
references there in). The subresultant techinques were introduced by Collins [H|, in order to 
reduce the growth of the coefficient in the signed polynomial sequences when pseudo-division is 
used. Gathen and Lucking m presented a unified approach to various definitions and algorithms 
of the subresultants, while El Kahoui m studied the subresultants in arbitrary commutative 
rings. For the Sturm-Habicht (or Sylvester-Habicht) sequences the reader may refer to the work 
of Gonzalez-Vega et al Emm 

In this paper we consider the Sturm-Habicht sequence of A and B, i.e StHa(A, B), which 
contains polynomials that are proportional to the polynomials in sPRS (A, B), i.e there exists an 
integer such that if we multiply a polynomial in StHa(A, B) we find the corresponding polynomial 
in sPRS (A, B). However Sturm-Habicht sequences achieve better bounds on the bit size of the 
coefficients and have good specialization properties, since they are defined through determinants. 
Moreover, they are obtained by sign modification of the subresultant sequence. 

Let Mj be the matrix which has as rows the coefficient vectors of the polynomials 

AX q ~’, AX q-2 ' -i ,..., AX, A, B, BX,..., BX p - 2 ^, BX”- 1 > 

with respect to the monomial basis X p+q_1_ l,X p+q_2_ l,... ,X, 1. The dimension of Mj is (p + 
q -1 -2j) x (p + q - 1 - j). 

For l = 0,..., p + q — 1 — j let Mj be the square matrix of dimension (p + q — 2j) x (p + q — 2)) 
obtained by taking the first p + q — 1 — 2j columns and the l-th column of Mj. 

Definition 2 The Sturm-Habicht sequence of A and B, is the sequence 

StHa(A, B) = (H p =H p (A,B),...,Ho=H 0 (A,B)) 

where H p = A,H p _i =B and Hj = Xu=o ^ e t (Mj )X l . 

The sequence of principal Sturm-Habicht coefficients 


(hp =hp(A,B),...,ho(A,B)) 


is defined as h. p = 1 and h.j = coef f j (Hj) is the coefficient ofX* in the polynomial Hj for 0 < i < p. 
When h.j = 0 for some j then the sequence is called defective, otherwise non-defective. 


If the StHa(A,B) is non-defective then it coincides up to sign with the classical subresultant 
sequence. The sign of proportionality is (—l) lp ,U 2 1 However, in the defective case, one has 
better control on the bit size of the coefficients in the sequence. 

One important property [2 is that the polynomial Ho(A, B), modulo its sign, is the resultant of 
A and B. Moreover the greatest common divisor of A and B is obtained as a by-product, together 
with the following equivalence: 


H k (A, B) = gcd(A, B) <H> 


ho (A, B) = ■ • • = Hk-i (A, B) = 0 
hk(A,B) ^0 


There are various algorithms that compute all polynomials in StHa(A, B), e.g EmUEZlEniEEl 
Most algorithms exploit the special structure of the cofficients of the polynomials that appear in 
the sequence and manage to reduce their bit size. All of these algorithms have more or less similar 
arithmetic and bit complexity. 
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Theorem 3 □ 03 77 Tfll There is an algorithm that computes StHa(A,B) in 0s(pq M (px)), 
or 0B(p 2 qx). Moreover, £(Hj(A,B)) = 0(p x). 

The complexity for the computation of the whole StHa(A, B) is optimal up to constants and some 
logarithmic factors. Notice that there are fl(q) polynomials in the sequence which have degree 
O(p), hence the total number of all the coefficients appeared in the sequence is fl(pq). This 
observation allows us to argue that the arithmetic complexity O(pq) achieved by the algorithms 
is optimal and since the bit size of the coefficients is O(px) one can also, trivially, deduce the 
optimality of the bit complexity. However, there are cases where only one polynomial in the 
sequence is needed (e.g some polynomial in the middle of the sequence or the gcd, or the resultant 
etc) or we do not need the actual sequence but the evaluation of it over a number. In these cases a 
faster algorithm exists, which represents the sequence implicitly by the quotient boot, and is based 
on a divide-and-conquer strategy and on the idea of the HALF-GCD algorithm. The idea of this 
“implicit” representation is not new. The reader may refer to Strassen EH, where the optimality 
of this evaluation scheme is proven basen on the work of Knuth and Shonhage and to Yap £23 
for the analysis of the HALF-GCD algorithm. Schwartz and Sharir 2E] mentioned the benefits 
of this approach for computations with real algebraic numbers and also Davenport El exploited 
the advantages of this approach for real root isolation. Lickteig and Roy m and independently 
Reischert El formulated this approach for Sturm-Habicht sequences. 

Theorem 4 H IM E3 fTp/ The quotient hoot, the resultant and the gcd of A and B. can be 
computed in 0n(q lg qM (px)) or 0b (p q x). 

Actually for the computation of gcd (A, B) various algorithms exist with complexity 0 b (pqx) (c.f 

EH El). 

Let the quotient boot that corresponds to StHa(A, B), be StHaQ(A, B) = (Qo, Qi,..., Qk-i. H k ). 
The number of coefficients in StHaQ(A, B) is 0(q) and their bit size is 0(p t) (c.f H ESI). The 
evaluation of the sequence on a number can be recovered from the quotient boot starting from 

H k . 

Theorem 5 There is an algorithm that computes the evaluation o/StHa(A,B) over a 

number a, where a £ Q U {± 00 } and has bit size at most cr in 0B(qlgqM (max (px, qcr))) or in 
0s(qM (max (px, qcr))) ?/StHaQ(A,B) is already computed. 

In both cases the complexity is 0b (q max (px, qcr)). 

Remark 6 In many cases (e.g real root isolation, comparison of algebraic numbers etc) we need 
the evaluation o/StHa(A,A ) over a rational number of bit size 0(pT). If we perform the evalua¬ 
tion by Horner’s rule then since there are Il(p 2 ) coefficients in the sequence and we must midtiply 
numbers of bit size 0(px) and 0(p 2 x), the overall complexity is 0b(p j M (px)). 

However, when we compute the complete StHa(A, A ) in 0b(p 2 M (px)) (Th.\E>, the quotient 
boot is computed implicitly J23DJ. Thus, we can use the quotient boot in order to perform the 
evaluation even if we have already computed all the polynomials in the Sturm-Habicht sequence. 

Remark 7 Notice also that the computation should be started by the last element of the quotient 
boot so as to avoid the costly computation of two polynomial evaluations using the Horner’s scheme. 

Even though this approach is optimal, we have to mention that involves big constants in its 
complexity, thus it is not efficient in practice if the length of the sequence is not sufficient big and 
special techniques should be used for its implementation in order to avoid costly operations with 
rational numbers. 

So, as it is always the case with optimal algebraic algorithms, the implementation is far from 
a trivial task. 

We can also use Sturm-Habicht sequences in order to compute the square-free part of A, i.e 
A re d = A/ gcd (A, A ), where A is the derivative of A. 
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Theorem 8 n Algorithm 10.17, p. 326] The square-free part of A, i.e. A red , can be computed 
from StHa(A,A), in Ob (p lgpM (px)) or (D B (p 2r r). Moreover, C(A red ) = 0(p + r). 

Remark 9 There is a normalization step m at the last element o/StHa(A, A ) so as to compute 
A re a which achieves the good bound for C (A red ) . Notice that if we rely on Mignotte's bound lw 
that applies to all exact polynomial divisors of A or on the subresultant algorithm then the bit 
size of A red is 0(p t). 

Let W(a,b)( 3 ) denote the number of modified sign changes of the evaluation of StHa(A,B) 
over a. Notice that W (AjB )(a) does not refer to the usual counting of sign varations, since special 
care should be taken for the presence of consecutive zeros EIE21I23|. 

Theorem 10 Let A, B € Z[X] be relatively prime polynomials, where is A square-free 

and A is the derivative of A. If a < b are both non-roots of A and y ranges over the roots of A 
in (a, b), then 


W (AiB )([a,b]) := W (A , B )(a) -W (A , B )(b) = sign(A'(y)B(y)). 

y 

Corollary 11 mu' If B = A then StHa(A, A ) is the Sturm sequence and Th. 1771 counts the 
number of real roots of A in (a, b). 

Remark 12 Actually Th. E3 computes the Cauchy index of the rational function B/A m 


3 Real algebraic numbers 

3.1 Univariate real root isolation 

Let f = Li-o a iN x € Z[X], with deg(f) = d and C (f) = x and let f re d be its square free part. We 
want to isolate the real roots of f, i.e to compute intervals with rational endpoints that contain 
one and only one root of f, as well as the multiplicity of every real root. 

Various algorithms exist for polynomial real root isolation, but most of them focus on square- 
free polynomials. The interested reader may refer to the algorithm of Collins and Loos 0, where 
the real roots of the derivative are used in order to isolate the real roots of the polynomial, 
with complexity O b ( d 9 + d 6 x 3 ), to the work of Akritas [T] for an algorithm based on continued 
fractions, or to the work of Rouillier and Zimmermann ':Tol (and references therein) where a unified 
approach with optimal memory management is presented for various algorithms that depend on 
Descartes’ rule of sign. The complexity of all the algorithms is no better than 0 B (d 6 x 2 ). Moreover 
Eigenwillig et al m. recently presented an algorithm for polynomials with bit stream coefficients 
which is based on Descartes’ rule of sign and the properties of the Bernstein basis. The complexity 
of their randomized algorithm is Ob (d 4 (log (sep) + x) 2 ), where sep is the separation bound (see 
Rem. E3>. that is 0[ 2 dT ), hence the complexity is 0 B [ d^x 2 ). 

However, we have to mention that the stated references are only the tip of the iceberg of the 
existing bibliography. 

If we restrict ourselves to real root isolation using Sturm (or Sturm-Habicht) sequences the 
first complete complexity analysis is probably due to Collins and Loos 0, that state a complexity 
of C>B(d 7 x 3 ). Davenport El improves this bound to 0 B (d 4 x 2 ) but does not present a formal 
proof that this bound holds for non square-free polynomials (Th. |SJ) and does not compute the 
multiplicities of the roots. Also Schwartz and Sharir m implicitly state this bound, but without 
a proof. Recently Du et.al m prove this bound for non square-free polynomials by giving a 
ingenious amortized-like argument for the number of subdivisions that must be performed. 

In this section we will prove that we can isolate the real roots of a non square-free polynomial 
in G B (d 4 x 2 ) and that in the same time we can also compute the multiplicities of the real roots. 


4 


Remark 13 We are not mentioning neither the algorithm of Mourrain et al that is roughly 
speaking based on a combination of Descartes’ rule and on the properties of Bernstein basis, nor 
its improvements I M since currently we are working with B. Mourrain on obtaining a complexity 
bound for this algorithm similar to the bound of the algorithm that uses Sturm-Habicht sequences. 

Algorithm 1 (Real Root Isolation using Sturm-Habicht Sequences) 

Input: A polynomial f € Z[X], with degf = d and C (f) = t. 

Output: A list of intervals with rational endpoints, which contain one and only one root oft 
and the multiplicitly over every real root. 

1. Compute the square free part oft, i.e f Te d 

2. Compute StHa(f Te d) 

3. Compute an interval Io = (—B,B) with rational endpoints that contains all the real roots. 
Initialize a queue Q with Io- 

4- While Q is not empty do 

Pop an interval I from Q and compute using Cor. m the number of roots in I. 

If I contains no real roots, discard I. 

If I contains one real root, output I. 

If 1 contains more than one real root split it to II and Ir and push them to Q. 

5. Determine the multiplicities of the real roots, using the square-free factorization of f. 

Remark 14 For a detailed description of Steps 2-4, that exploits the implementation details, the 
reader may refer to Davenport et al m or the more recent work of Du et al m ■ Notice that special 
care should be taken for the case that the middle of a tested interval is a root of the polynomial, 
which is a non trivial implementation issue. 


3.1.1 Complexity analysis of real root isolation 

Step 1 The computation of f Te d can be done in C>B(d 2 T) (Th. 0). Notice that C (f Te d) = O (d + T). 
We assume that d = 0(r), thus C (f Te d) = C?(t). 


Step 2 We do no need the complete sequence (Remark 0, we only need the quotient boot, thus 
this computation can be done in Op, (d 2 T) (Th.QJ. However, we may also assume that the 
complete sequence is computed, with complexity Oe(d 3 T) (Th.0), since this step is not the 
bottleneck of the algorithm. 


max 


Step 3 The Cauchy bound [U EH EH states that if a is a real root of f then |a| < B = 1 + 
— , a , f ‘ 2 ,..., 77 s - ). Various upper bounds are known for the absolute value of 

In a n &n J 

the real roots (c.f DEI 1321101). However, asymptotically the bit size of all the bounds is 
the same, i.e B < 2 T . 


Step 4 We count the number of real roots using Cor.^] by evaluating StHa(f re d) over rational 
numbers of bit size at most O (dx) CRemark l20H . The cost of every such evaluation is O-q (d 3 t) 
(Th.0. Since the number of subdivisions that we must perform in order to isolate all the 
real roots is 0(dx + dlgd) (Prop. EJ, the overall complexity of this step is C>B(d 4 T 2 ). 
Notice that the complexity of this step dominates the complexities of all the other steps. 
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Step 5 In order to compute the multiplicities we compute the square-free factorization, i.e a sequence 
of square-free coprime polynomials (gi, g 2 ,..., g m ) with f = gi g\ ■ ■ ■ g™ and g m k 1. The 
algorithm of Yun m computes the square free factorization in (Db (d 2r r). To be more specific 
the cost is twice the cost of the computation of StHa(f, f ) m- 

At every isolating interval, one and only one g, must have opposite signs at its endpoints, 
since gj are square free and pairwise coprime. If gj changes sign at an interval then the 
multiplicity of the real root that the interval contains is j. 

We can evaluate each gi, 1 < i < m, at all the isolating points simultaneously in Oe(d 3 T) 
H3E|. To be more precise we can perform an evaluation of a polynomial at d numbers at 
the cost of evaluating the polynomial over one number using Horner’s scheme. Since m is 
at most d, the overall cost is C>B(d 4 T). 

This bound is quite pessimistic due to the reason that it is not possible for the square free 
factorization of f to contain O(d) polynomials of degree 0( d) since this would lead to bound 
of 0(d 2 ) for the degree of f. 

Thus we may assume that either m is a constant, or that the degrees of gj’s are bounded 
by a constant. Hence the complexity of this step is 0^(d ir x). However there is no need for 
a detailed study of the complexity of this step since both mentioned complexities do not 
dominate the complexity of the overall algorithm. 

To complete the proof of the complexity of the algorithm we must prove that the num¬ 
ber of subdivisions is 0 [dT + dlgd). Recall [2] that the Mahler’s measure, of f is M(f) = 
K\nt 1 max{l, |yi|}, where q p is the leading coefficient and Yi are all the roots of f. We know 
that Mlf) < 2 T ^^^ + 7 0. 

Remark 15 Since the product of the absolute values of roots greater than 1 of f re d can not be 
greater than the one off, we have the following inequality m 

At (fred) < M[f) < 2 T V / d~+T 

For the minimum distance between consecutive real roots of a square free polynomial the 
Davenport-Mahler bound is known E!. The conditions of this bound where generalized by Du et 
al (13) . However, using Remark ll5l we can provide a similar bound for non square free polynomials 

El- 

Theorem 16 (Davenport-Mahler bound revisited) Let ai < ot 2 < • • • < <Xk < (Xk+i be the 
k + 1 distinct real roots of f , which is not necessarily square free, with k > 1. Then 

k nr 

M(f) > n l«i - Oi+11 > -M(f)- d+1 d-T(Ai) k 

t=i d 

Proof. If f is square free then the bounds hold El 

If it is not square free let n < d be the degree of f T ed- Notice that k + 1 < n < d. Since the 
bound holds for f re d, we have 


M( fred) 

> 

nL 

kt 

- Ot+1 1 

> 

Ml fred)- n+1 n-5(^) k 


M{ f) 

> 

nti 

kt 

- 0ti+1 

> 

7W(f)- n+1 n'T(^) k 

(Rem. ESJ 

M{ f) 

> 

nL 

kt 

- OC i+ 1 1 

> 

Ad(f)- d+1 d-T(^) k 

(n < d) 


Thus the theorem holds. □ 

Proposition 17 The number of subdivisions that we need to perform in order to isolate the real 
roots of fred using Sturm-Habicht sequences is at most (D[ dT+ dlgd). 
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Proof. The proof follows Ei and uses the result of th. m 

Let oci be the real roots of f re d in increasing order. We know that the roots are contained in 
an interval I = (—B, B) (we can compute such an interval using for example the bound of Cauchy, 
i.e B < 2 T [2 [7 I40|h If f re d has k + 1 real roots then we need to find k separation points. For 
two consecutive real roots oq and oq+i a separation point for them is of magnitude j\a\ — ai + i| 
and we can compute it with at most [lg a j ] subdivisions, using binary search. 

Let 5(1) denote the total number of subdivisions that we need to perform in order to compute 
all the isolating points of f Te d in I- Then 


sm <Y_ 


i = l 


!gi 


2B 


I at - oci+il 


k 

< k + klg(2B) -^lg|at - a i+1 | 

i=1 


where the additional k in the last inequality represents the k possible roundings. 

5(1) < k + klg (2B) + (d — 1) lg M(i) + j lgd — klg\/3 + klgd (Th. Uhl) 

< k + klg (2B) + (d — 1) lg (2 T Vd+T) + -j lgd — klg\/3 + klg d (Rem. H5ll 

< k + k(x + 1) + (d — 1 )t + 2dlg d + klg d (B = 2 T ) 


Furthermore, if k < d — 1 =1> k < d then 5(1) = 0(dr + dig d). 


□ 


Remark 18 Du et al m obtained a similar bound using a charging scheme (amortized analysis) 
for each subdivision that may provide better constants and that can also be used for complex root 
isolation. 

Under the hypothesis that d = O(t) (to simplify notation), the previous analysis leads to the 
following theorem 

Theorem 19 Let f £ Z[X], with deg(f) = d and 5(f) = t, not necessarily square-free. We 
can isolate the real roots of f and determine their multiplicities using Sturm-Habicht sequences 
in OB(d 4 T 2 ). Moreover, the numerator and the denominator of the endpoints of the isolating 
intervals have bit size bounded by 0( dr). 

Remark 20 The bit size of the endpoints of the isolating intervals is bounded by the bit size of the 
separation bound off, i.e the minimum distance between two consecutive real roots. It is known 
PI [23 that sep(f) = sep(f re d) > d _i ^ (d + 1 j 1 ? 1 2 T(1 ~ d) , thus lg(sep(f)) = lg(sep(f re d)) = 
0[ dr). 

If a separation bound only on the real roots is needed then a bound due to Rump may be 
used, i.e sep(f) > d \/8^1 +max|ai| d ^ , which is a bit sharper but asymptotically has the 
same bit size. 

Remark 21 An + lg d) 2 ) bound on the complexity for univariate real root isolation can 

be obtained if we adopt the algorithm of Mourrain et al JM I. which is based on Bernstein basis and 
seems to have the best complexity in practice. The interested reader may refer to HEF for more 
details. 

3.2 Representation of real algebraic numbers 

The real algebraic numbers, i.e. those real numbers that satisfy a polynomial equation with integer 
coefficients, form a real closed field denoted by K a ig = Q- From all integer polynomials that have 
an algebraic number oc as root, the one with the minimum degree is called minimal polynomial. 
The minimal polynomial is unique, primitive and irreducible HUES In our approach, since we 
use Sturm-Habicht sequences, it suffices to deal with algebraic numbers, as roots of any square-free 
polynomial and not as roots of their minimal ones. 

In order to represent a real algebraic number we chose the isolating interval representation. 
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Definition 22 The isolating-interval representation of real algebraic number a £ R Q i g is a = 
(P(X),I), where P(X) £ Z[X] is square-free and P(a) = 0, I = [a,b], a,b,G Q and P has no other 
root in I. 

Remark 23 Using the results of Section lM.il we conclude that we can represent all the roots of a 
polynomial f, with deg(f) = d and C (f) = x, using isolating interval representation in C B (d 4 x 2 ) 
and that the endpoints of the isolating intervals have bit size 0( dx). 

This will be the case for all the real algebraic numbers that we will consider for the rest of the 
paper, except stated otherwise. 

3.3 Comparison and sign evaluation 

We can use Sturm-Habicht sequences in order to find the sign of a univariate polynomial, evaluated 
over a real algebraic number and to compare two algebraic numbers (cf. |17l for degree < 4, where 
it is proven that these operations can be performed in 0( 1)), or 0 B (x)). 

Lemma 24 E3EIEI Let Q(X) £ Z[X], where deg(Q) = d and C (Q) = T, and a real algebraic 
number a = (P, [a, b]) . We can compute sign(Q( a)) mO B (d 3 x). 

Proof. By th.EJJ sign(Q(a)) = signfW^Q [a, b] • P (a)). Thus we need to perforce two evaluations 
of StHa(P, Q) over the endpoints of the isolating interval of a. The complexity of each is Ob (d°x) 
(Th. 0 and Rem. EH- which is also the complexity of the operation. □ 


Lemma 25 030 ' We can compare two real algebraic numbers in isolating interval representa¬ 
tion in C B (d 3 x). 

Proof. Let two algebraic numbers yi = (Pi(x),Ii) and yi — (PiM.O where Ii = [ai,bi], 
I2 = [a2, b2]- Let J = Ii n I2. When J = 0 , or only one of yi and yx belong to J, we can easily 
order the 2 algebraic numbers. If yi ,y2 £ J, then yi > yi 44 P2(yi) ■ P 2 (y2) > 0. 

We obtain the sign of P 2 (y2), using Lem. [221 thus the complexity of comparison is 0 B (d 3r r). 

□ 


3.4 Simultaneous inequalities 

Let P, Ai,..., A n ,, Bi,..., B n2 , Ci,..., C n3 £ Z[X], with degree bounded by d and coefficient bit 
size bounded by x. We wish to compute the number and the real roots, y, of P such that Ar(y) > 0, 
Bj (y) < 0 and Cic(y) = 0 and 1 < i < ni, 1 < j < n. 2 ,1 < k < n. 3 . Let n = ni + 112 + 1 x 3 . 

Corollary 26 There is an algorithm that solves the problem of simultaneous inequalities time 
Or (max{n d 4 T, d 4 x 2 }) = Or (d 4 x max{n, t}) . 

Proof. First we compute the isolating interval representation of all the real roots of P in Or (d 4 x 2 ) 
(Th.GH). There are at most d. For every real root y of P, for every polynomial At, Bj, Ck we 
compute the sign (At (y)), sign (Bj (y) ) and sign (Ck(y)). 

Sign determination costs Or[ d 3 x) (Lem. m and in the worst case we must compute n of 
them. Thus the overall cost is C B (max{nd 4 x, d 4 x 2 }). □ 


Remark 27 Ben-Or, Kozen and Reif considered the problem of simultaneous inequalities. 
The interested reader may also refer to the work of Canny for a variant that is faster in the 
univariate case (it saves a factor) which has arithmetic complexity 0 (u(mdlg mlg 2 d + m 2 " 376 )), 
where m is the number of real roots o/P. 


Coste and Roy m introduced Thom's encoding for the real roots of a polynomial and the prob¬ 
lem of simultaneous inequalities in this encoding. Their approach is purely symbolic and works 
over arbitrary real closed fields. They also derive m polynomial bounds for the arithmetic com¬ 
plexity of the problem. This is not the case for our approach, which is pseudo-polynomial as all 
the approaches that depend on separation bounds, since for the real root isolation the number of 
subdivions that we must perform depends on the bit size of the coefficients. They state a com¬ 
plexity of Ob[t a 8 ), using fast multiplication algorithms but not fast computations and evaluation 
of polynomial sequences, where m is an integer bigger that u, d, t. In this notation our bound is 
C>B(vn 6 ). 

We are planning to perform a more thorough analysis of algorithms that use Thom’s encoding by 
embedding the results of fast multiplication and fast computation and evaluation of Sturm-Habicht 
sequences, as well as an efficient implementation. 

Remark 28 In the book of Basu, Pollack and Roy m an algorithm for the problem is pre¬ 
sented when the real algebraic numbers are in isolating interval representation, with complexity 
O B(n.d 6 T 2 ). However their analysis does not assume fast multiplication algorithms and they de¬ 
termine the sign of a polynomial over an algebraic number by repeated refinements of the isolating 
intervals. 

3.5 Computation in an extension field 

If a = (A, I), then we can perform computations in Q(<x). If deg(A) = d then Q(a) is a vector 
space with basis elements (1, a,..., a d_1 ) m Thus we can represent every element |3 £ Q(a) as 
polynomial of degree at most d — 1 with integer coefficients (this is wlog since we can always clear 
denominators). 

The four elementary operations in Q(<x) are the usual operations with polynomials, however 
special care should be taken m if A is not the minimal polynomial of a, which is usually the 
case. The complexity of the four operations are the complexity of the corresponding polynomial 
operations. However we can improve a lot the practical complexity if we represent the elements of 
Q(a), which are univariate polynomials, in the Horner’s basis (Jj. The latter approach allow us to 
argue that the complexity of the four operations is quasi-linear, hence optimal up to logarithmic 
factors. 

The most difficult elementary operation is the determination of the sign of 3 £ Q(a), which 
also corresponds to the comparison of two numbers that belong to the same extension field. Since 
P is represented by a polynomial P, it suffices to compute sign(P(a)). This can be done easily 
using the results of Lem. 1771 

Remark 29 Currently we are implementing an Algebraic_extension class in SYNAPS VT(t . This 
class allows us to perform real root isolation and other elementary operations with polynomials 
that belong to Q(a)[X]. It is not clear at all, at least from a practical point of view, whether 
the computations of the signs of various quantities that are needed should be performed using 
Sturm-Habicht sequences or repeated refinements of the isolating intervals. We believe that a 
combination of both approaches will eventually lead to an optimal, from an implementation point 
of view, scheme. 

Remark 30 Computational, is very costly to perform operations with two numbers that belong 
to two different extension fields, or to change the extension field of a number. Currently only the 
package of Rioboo in Axiom HF can perform non trivial operations with real algebraic numbers 
that belong to different extension fields (actually to towers of extension fields). 

However given an extension field Q(<x) we can easily compute the representations of the Q(—a), 
Q(^), Q(a ± a), Q(a of) and Q(-f)? a £ Q, where by representation we mean the isolating interval 
representation of the real algebraic number that defines the extension field. 

We will present the full details of algorithmic, complexity and implementation issues for com¬ 
putation in an extension filed(s) in a future report since this is work in progress. 
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Remark 31 Recently we have implemented in MAPLE a prototype library, that besides real root 
isolation and sign evaluations also provides the four operations between real algebraic numbers that 
belong to different extension fields as well as fractional powers of real algebraic numbers. In the 
near future we will make the library freely available. 

Remark 32 As for C++ implementations, with freely available code, that can perform the basic 
operations with real algebraic numbers, i.e real root isolation and sign evaluations, the reader 
may refer to SYNAPS TW . where also the bivariate problems of the next section are treated, or 
to the library of Guibas et al \2f\ , especially optimized for kinetic data structures, or to NiX the 
polynomial library of EXACUS & which is part of a bigger library for non linear computational 
geometry. 

4 Computations with 2 real algebraic numbers 

In this section we will present some algorithmic and complexity results concering computations 
with two real algebraic numbers and real solving of bivariate polynomial systems. 

The algorithms, the implementation details as well as experimental results for bivariate poly¬ 
nomial system solving were presented in GEi. In this work we present the complexity analysis. 

4.1 Sturm-Habicht sequences for bivariate polynomials 

Theorem 33 Ifffl Let F, G £ (Z[Yi,..., Y t ])[X], where deg x (F) = p > q = deg x (G) and C (F) = 
C (G) = t. Moreover, let deg Y (F) < &i and deg Y . (G) < &*, where I < i < t. We can compute the 
quotient boot of F and G, the resultant and the gcd in 

O B (qlgqM((p + q) t+ 1 6 1 •••5 t T)) 

Remark 34 Th. [23 generalizes Th. It is easy to deduce this complexity from the univariate 
case if we transform the multivariate polynomials to univariate ones, using iterated applications 
of univariate “binary segmentation” and Kronecker’s map :3 Ef- In an expanded version of 
this paper we will present complexity results for multivariate Sturm-Habicht sequences based on a 
more sophisticated binary segmentation that can save some logarithmic factors and that has 
a simpler encoding and decoding algorithm, that involves only shifts and additions. 

Remark 35 If F and G are bivariate polynomials the complexity of computing the quotient boot, 
the resultant and the gcd is O b (q lg qM (r(p + q) 2 &i)) = O b (qp 2 5i t). The bit size of the gcd as 
well as the bit size of the polynomials in the quotient boot is bounded by 0 (pT) 0E3EI- 

The fact that Sturm-Habicht sequences are amenable to any specialization of the coefficients 
EH 1221 finds application when we are computing with multivariate (and in our case bivariate) 
polynomials. 

Let F and G be two polynomials with parametric coefficients, such that their degree does not 
change after any specialization in the parameters. The computation of their Sturm-Habicht se¬ 
quence before specialization of their coefficients, guarantees that the seuenece is valid under every 
specialization. We use this property so as to compute such a sequence for bivariate polynomi¬ 
als, regarding them either as polynomials in (Z[X])[Y] or in (Z[Y])[X], Remember that the last 
polynomial in the sequence is the resultant with respect to X or Y, respectively. 

The following theorem will be very useful 

Theorem 36 HD2JQSI Let f, g square-free and coprime polynomials, such that Cf and C g are 
in generic position. If 

Hj(X,Y) =StH aj (f,g) =h j (X)Y i +h jJ _ 1 (X)Y i - 1 + ••■+h,, 0 (X) 
then if C = (oc, (3) £ Cf (~l C g then there exists k, such that 

h 0 (oc) = ••• = h k - 1 (a) = 0 , h k (a)^ 0 , 

k h k (a) 
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4.2 Bivariate sign evaluation 

The previous tools suffice to compute the sign of a bivariate polynomial function evaluated over 
two algebraic numbers. Consider F G Z[X,Y] and oc = (A(x),Ii) and |3 = (B(X), I 2 ) where 
Ii = [ai,bi], I 2 = [ci 2 ,b 2 ]. We wish to compute the sign of F(a, |3). 

Theorem 37 (Bivariate sign_at) Let F G Z[X,Y] with deg x (F) = deg Y (F) = di and £ (F) = x 
and two real algebraic numbers a = (A, I a ), (3 = (B, Ip) where A, B G Z[X], deg(A) = deg(B) = d 2 , 
£ (A) = £ (B) = t and I<x, Ip G Q 2 with bit sizes bounded by 0(6.2 x). 

We can compute the sign of F evaluated over oc and (3, i.e sign(F(a, |3)) in Ob (df d| t) time, 
assuming di < d 2 . 

Proof. In order to compute the sign of F(ct, (3), we consider F as a univariate polynomial in X, 
i.e F G (Z[Y])[X] and we try to compute its sign when we evaluate it over oc, as we did in the 
univariate case (Cor.^J). 

We compute the Sturm-Habicht quotient boot of A and F with respect to X in (^(dld^x), as 
follows: 

We transform F to a univariate polynomial in Z[X], using binary segmentation induced by the 
map 4> : Y 1—> 2 C d2T [231231101, where c is a suitable constant. Now £ (F) = 0 [ di d 2 x). Thus, the 
computation of StHaQ(A,F) takes d?B(d?d 2 x) time, using Th. 0 since p = di , q = d 2 and the 
bit size is 0(6i d 2 x) and dominates the cost of binary segmentation. 

Notice that there are 0( di) polynomials, in Z[X], in the sequence. 

We evaluate the sequence on the left endpoint of I a , which has bit size 0 (d 2 x). The evaluation 
of the sequence can be done in C>B(d 2 d 2 x), using Th. |3 where p = d 2 , q = di, cr = 0(6ix) and 
the bit size is 0 (di d 2 x). 

By applying the inverse transformation of 4> in the evaluated sequence we get 0( di) univariate 
polynomials with respect to Y, with degrees 0 (did 2 ) and bit size C>(did 2 x). For every such 
polynomial we compute its sign when we evaluate it over |3. This can be done in Ob (dfd^x) (Th.[SJ), 
which dominates the cost of the inverse transformation. Since there are 0( di) polynomials, this 
step can be done in Ob( d^d^x). 

We do the same for the right endpoint of I a . This suffices to compute sign(F(a, (3)) using 
Th.^ilor Cor. 1241 

Thus the overall complexity of the operation is Ob (df dfx). □ 


Remark 38 The complexity of computing StHaQ(A,F) can also be derived from Rem.. 17751 where 

p = d 2 , q = 61 = di. 

Remark 39 We can extend this approach to polynomials with an arbitrary number of variables, 
similar to However the usage of Sturm-Habicht sequences, instead of generalized Sturm 

sequences, improves both the theoretical m and the practical complexity {HEBEI- 

4.3 Two variants of bivariate real solving 

In what follows we will present two variants of bivariate polynomial real solving. We assume that 
the polynomials are square-free. However we can easily drop this assumption. 

4.3.1 Modified RUR 

We use Tlil.'idl following |231 GU, so as to compute the solution of bivariate polynomial systems. 
We consider polynomials f, g G Q[X, Y], such that Cf,C g are in generic position and we compute 
the resultant of f, g with respect to Y, which is a polynomial in X. The real solutions of the 
polynomial correspond to the x— coordinates of the solution of the system. Then, using Th. ESI 
we lift these solutions in order to determine the y— coordinates, as a rational univariate function 
evaluated over an algebraic number. Even though the previous approach is straightforward, it has 
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one main disadvantage. The ^-coordinates are computed implicitly. If this is all that we want 
then this is not a problem. However in most cases we want to further manipulate the solutions 
of the system, i.e. to compare two y—coordinates or to count the number of branches of each 
curve above or below this ordinate. Of course we can always find the minimal polynomial of these 
algebraic numbers, but this is quite expensive. Thus we chose an alternatively way. 

We compute the resultant, using the Sturm-Habicht Sequence, both with respect to Y and X, R x 
and Ry respectively. We compute the isolating interval representation of the real roots of R x and R u 
(Th.EHJ Let cxi < • ■ ■ < <Xk and (3 1 < • • • < (3 r be the real roots of R x and R y , respectively. For the 
real roots of R y we compute rational intermediate points, q o < (3 1 < qi < • • • < qi_i < (3 1 < q r 
where qj £ Q, 0 < j < l. We can easily compute the intermediate points, since the algebraic 
numbers are in isolating interval representation. 

For every root oq, 1 < i < l, using Th. m we compute a rational univariate representation of 
the corresponding y-coordinate, which is without loss of generality, of the form yi = 4 . Since 

have already computed the real solutions of R y , it suffices to determine to which (3j, yi equals to, 
that is to find an index j such that 


qi < 


A(aQ 

B(ai) 


< qj+i 


or, if we assume that A(at) > 0, this can be checked using Lemma EH then 


qjA(ai) < B(oq) < q j+1 A(a t ) 

Actually what we really want is to determine the sign of univariate polynomials of the form 
U(X) = qjA(X) — B(X) evaluated over the real algebraic numbers that are solutions of R x = 0. 
This can be done using Lemma, PHI 

Remark 40 We computed the complexity of the modified RUR algorithm is 0B(d lo T 2 ), since the 
complexity of the algorithm is dominated by the real solving of the univariate solutions. However 
we will present the complexity analysis in an extended version of the paper, since we are currently 
working on the complexity of transforming two algebraic curves to generic position and on the 
complexity of transforming the real solutions back to the original coordinate system if a random 
shear is performed as well as on the complexity of the topology of real algebraic curves (see also 
Rem. \fl\) . 

Remark 41 (Topology of an algebraic curve in 2D) As stated by G. Vega and El Kahoui 
m the bottleneck of the algorithm for the topology of a plane real algebraic curve is the computation 
of the Thom’s codes of the real roots of a univariate polynomial. In our notation this is O b (d 1 4 t 2 ). 

The complexity results that we presented for univariate real root isolation fTh. il hi) can improve 
their bound. We claim that the bound for the topology algorithm is O^(d^ 0 i: 2 ), induced by the 
modified RUR algorithm, even when the curves are not in generic position. However, for this very 
crucial problem in non-linear computational geometry a more detailed analysis is needed so as to 
give light to various details of the algorithm in order to present an efficient implementation for 
real-life applications. 

4.3.2 Naive real solving 

In many cases, especially in non-linear computational geometry or when we want to solve a system 
of inequalities, it is very difficult to deal with bivariate polynomials that are not in generic position. 

Even though the assumption of generic position is without loss of generality since we can 
apply a transformation of the form (X,Y) i—» (X + aY, Y), where a is either a random number 
or a number computed deterministically EDI before the execution of the algorithm, or we can 
detect non-generic position during the execution Eau, then apply a transformation of the form 
(X,Y) i—» (X + Y, Y) and start the algorithm recursively. However neither approach is an easy 
computational task. Moreover if we need the solution in the original coordinate system then we 
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must perform the inverve transformation, and this is a very difficult task, especially when we want 
to avoid refinements that can go up to seperation bounds. 

In order to overcome such barriers we suggest one (naive) variant for bivariate polynomial 
system solving. The intersting reader may refer to m for a complete descritption of the algorithm 
and for preliminary experimental results. 

Theorem 42 (Naive solve) Let F, G G Z[X, Y] be square-free polynomials with total degree bounded 
by d and (coefficient) bit size bounded by t. We can compute the real solutions of the system 
F = G = 0 in e> B (d 14 T), under the assumption that t = 0 B (d 4 ). 

Proof. We compute the resultant of F and G with respect to Y, i.e Rx in time 0 B (d 4 x) (Rem. l35l) . 
Note that deg(R x ) = deg(Ry) = 0[ d 2 ) and C (Rx) = £ (Ry) = 0{ dx). 

We compute the real algebraic numbers that are roots of Rx, in isolating interval representation 
in 0 B (d lo x 2 ) time (Th. ITTH) . Their representation involves the square-free part of Rx, also of 
coefficient bit-size (9(dx) (Th. [SJ), and an isolating interval, whose endpoints have bit-size 0( d 3 x) 
(Rem. l2()ll. There are at most 0( d 2 ) real roots, i.c deg(Rx). We do the same with respect to Y. 

We form all the possible pairs of real algebraic numbers and for every pair of algebraic numbers 
we test if both F and G are zero. There are G( d 4 ) such tests and each takes C> B (d 1 °x) time 
(Th. Ef 

Thus the overall complexity is 0 B (d 14 x). □ 


4.4 Bivariate simultaneous inequalities 

Let P, Q, Ai,..., A n ,, Bi,..., B n2 , Ci,..., C n3 € Z[X, Y], with degree bounded by d and coef¬ 
ficient bit size bounded by x. We wish to compute the number and the real roots, (a, 3), of 
the system P(X,Y) = Q(X,Y) = 0 such that At(a, |3) > 0, Bj(a, (3) < 0 and Ck(a, (3) = 0 and 
1 < i < ni, 1 < j < t \2 ,1 < k < 1 x 3 . Let n = ni + ri 2 + 1 x 3 . 

Corollary 43 There is an algorithm that solves the problem of simultaneous inequalities in time 
Ob (max{nd 12 T, d 14 x}) = 0 B (d 12 xmax{n, d 2 }). 

Proof. First we compute the isolating interval representation of all the real roots of the system 
P = Q = 0 in C> B (d 14 T) (Th. El- There are at most d 2 real roots of the system. For every real 
root (a, (3) of P, for every polynomial At, Bj, Ck we compute the sign At(a, (3), sign (Bt(a, (3)) 
and sign (Ci(a, |3)). 

Sign determination costs C> B (d 10 T) (Th. 13171 and in the worst case we must compute d 2 of 
them. Thus the overall cost is C B (max{nd 12 T, d 14 x}). □ 


5 Conclusions and future work 

We plan to apply our tools in computing the topology of algebraic curves in 2D and 3D, as well 
as the topology of surfaces in 3D. Another possible approach, to be implemented and compared 
at a practical ant theoritical level, includes the adoption and Thom’s encoding [annum Last 
but not least, we intend to use arithmetic filtering to handle cases that are far from degenerate, 
so as to improve the speed of our software for generic inputs. 
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